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1. 


INTRODUCTION 


The  ability  of  modern  instrumentation  to  collect 
( local  gravity  data  is  quickly  surpassing  our  capability  to 
process  and  interpret  this  data.  For  example , the  NASA  GEOS-3 
satellite  has  already  gathered  20  million  measurements  of 
geoid  height.  Storage,  analysis,  and  evaluation  of  this 
data  is  a formidable  task. 

Fourier  transforms  are  an  efficient  tool  for  analy- 
zing such  bulk  data  because  space-domain  convolutions  are  re- 
placed by  frequency-domain  multiplications.  In  principle, 
Fourier  transforms  can  be  used  for  almost  all  data  processing 
tasks:  upward  continuation,  computing  anomalies  from  undula- 

tions, smoothing,  least-squares  collocation,  etc.  However, 
Fourier  methods  are  based  on  the  flat-earth  approximation. 

The  accuracy  of  this  approximation  is  therefore  a key  issue. 
Flat  earth  theory  is  only  valid  on  a local  basis,  but  how 
large  is  a "local"  area?  Do  flat-earth  methods  always  lead 
to  long-wavelength  errors? 

The  purpose  of  this  paper  is  to  develop  improved  flat- 
earth  approximations  and  thereby  enhance  the  usefulness  of 
Fourier  methods  in  physical  geodesy.  These  improvements  are 
accomplished  by  treating  flat-earth  theory  as  the  leading  term 
of  an  asymptotic  expansion. 

SRE  ~ SFE  + £g2  + £ S3  + • • • (1-1) 


Here  gRtr  and  gpp  are  round-  and  flat-earth  gravity  models. 
The  remaining  terms  on  the  right  hand  side  are  first-  and 


higher-order  corrections  to  flat-earth  theory.  The  small 
parameter  z is  inversely  proportional  no  the  earth’s  radius 
and  goes  to  zero  in  the  limit  of  large  earth  radius.  As  will 
be  shown,  the  correction  terms  are  relatively  easy  to  calcu- 
late using  Fourier  transforms.  Therefore,  the  flat-earth 
approximation  can  be  improved  to  any  desired  accuracy.  In 
practice , the  first  correction  term  is  often  adequate  because 
the  parameter  £ is  small  (s<<l).  Loosely  speaking,  the  para- 
meter z is  small  whenever  the  gravity  signals  of  interest  are 
concentrated  in  wavelengths  that  are  short  relative  to  the 
radius  of  the  earth.  This  situation  usually  prevails  because 
90  percent  of  the  energy  (variance)  in  gravity  anomalies  is 
contained  in  wavelengths  shorter  than  3000  km  (Ref.  1). 

For  typical  problems  in  geodesy,  the  ("inner")  ex- 
pansion (1)  converges  only  in  the  vicinity  of  the  point  where 
the  flat-earth  approximation  has  been  applied.  However,  this 
nonuniformity  can  be  corrected  by  forming  an  "outer"  expansion 
that  is  valid  at  large  distances.  These  expansions  are  com- 
plementary, share  an  overlapping  region  of  convergence,  and 
can  be  combined  (matched)  to  produce  a unif ormly-valid  rep- 
resentation of  both  the  local  and  global  field.  In  this  paper, 
inner  and  outer  expansions  are  derived  and  appropriate  matching 
procedures  are  given. 

The  approach  is  based  on  the  theory  of  matched  asymp- 
totic expansions  (Refs.  2,  3).  This  theory  was  originally 
developed  for  boundary-layer  problems  in  fluid  mechanics  but 
has  found  applications  in  many  other  areas  of  applied  math- 
ematics. In  the  present  context,  flat-earth  geodesy  is  the 
first  term  of  an  inner  expansion  (Eq.  1-1). 

The  connection  between  physical  geodesy  and  boundary- 
layer  theory  is  not  accidental.  Gravity  anomalies  are  caused 
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chiefly  by  terrain  and  the  associated  isostatic  compensation 
at  the  Mohorovicic  discontinuity.  Most  of  the  energy  in  ter- 
rain is  concentrated  in  wavelengths  that  are  short  relative 
to  earth's  radius.  Also,  the  earth's  crust  is  relatively 
thin  compared  to  the  earth’s  radius.  Therefore,  a "gravity 
boundary  layer"  exists  around  the  earth  with  a thickness  of 
a few  hundred  kilometers.  Within  this  boundary  layer  the  flat- 
earth  approximation  and  Fourier  transforms  are  appropriate. 
Outside  the  boundary  layer,  the  curvature  of  the  earth  is  im- 
portant and  spherical-earth  formulas  are  needed.  The  theory 
of  matched  asymptotic  expansions  provides  a systematic  method 
of  matching  solutions  that  are  valid  inside  and  outside  the 
boundary  layer. 

In  the  next  two  sections  the  method  of  matched  ex- 
pansions is  illustrated  with  examples.  These  examples  pro- 
vide an  introduction  to  the  method  but  are  also  important 
in  their  own  right  because  they  give  the  appropriate  outer 
expansions  that  are  needed  in  later  sections  of  the 
paper . 


1.1  DETERMINISTIC  EXAMPLE  - SHALLOW  POINT  MASS 

The  simplest  illustration  of  the  flat-earth  approxi- 
mation is  a point  mass  buried  at  a depth  D that  is  shallow 

relative  to  the  radius  of  the  earth  R (D<<R  ).  The  distur- 

e e 

bance  potential  T is 

T(r,i|0  = y jjr2  + ( 1-e ) 2 - 2r(  1-e)  cos  dT!  _1/ 2 (1.1-1) 

where 

r » r’/R0  e = D/Re  (1.1-2) 

u = km/Rg  (1.1-3) 


Here  r'  and  'jj  are  the  radius  and  angle  from  the  center  of  the 
earth  (Figure  1.1-1).  Also,  k and  m are  Newton's  gravitational 
constant  and  the  mass  of  the  disturbing  body.  The  Legendre 
expansion  of  Eq.  1.1-1  is 


T(r,t|/) 


l 

n=0 


(cost; ) 


(1.1-4) 


R-3046S 


Figure  1.1-1  Notation  for  Buried  Point  Mass  Example 


Because  the  parameter  s is  small,  the  Legendre  series  Eq . 1.1-4 
converges  slowly  and  is  inaccurate  for  numerical  evaluations. 

A more  useful  series  is  obtained  by  expanding  in  powers  of  the 
parameter  z: 


(r,v)  ~ 


(r2^l-2r  cost; ) 


1-r  cos  'li 

9 

r J-l-2r  cost; 


(1.1-5) 


However,  this  series  fails  to  converge  in  the  vicinity  of  the 
point  mass  (r<l+e  and  as  can  be  seen  from  the  ratio 

test.  In  the  terminology  of  perturbation  theory,  the  asymp- 
totic expansion  (Eq.  1.1-5)  is  "aonuniform"  and  the  perturbation 
problem  under  consideration  is  "singular"  rather  than  "regular". 
Equation  1.1-5  is  referred  to  as  an  "outer"  expansion  hence  a 
subscript  (TQ)  is  used. 

The  standard  remedy  for  such  a nonuniform  expansion 
is  to  define  "inner"  variables  that  magnify  the  region  of 
nonuniformity.  This  leads  to  an  inner  expansion  that  is  valid 
in  the  vicinity  of  the  point  mass.  In  general,  the  appropriate 
magnification  is  dictated  by  the  characteristic  dimension  of 
the  nonuniformity.  In  this  example  the  altitude  (r-1)  and 
angle  both  need  to  be  stretched  by  the  factor  t.  Thus,  inner 
variables  (z,R)  are  defined  as 

z = (r-l)/s  and  R = ty/z  (1.1-6) 


The- inner  expansion  (Eq.  1.1-7)  is  obtained  by  substituting 
Eq.  1.1-6  into  Eq . 1.1-1 


Tt(z,R) 




C(l-z)2 


R2]1'2 


( 1-z ) R2 
(1+z)2  + R 


2 


Cl. 1-7) 


The  first  term  of  the  inner  expansion  is  the  familiar  flat- 
earth  approximation.  3y  computing  the  next  term  in  the  series 
(Eq.  1.1-7)  it  can  be  shown  that  the  series  diverges  for 
R>0(l//t),  or  equivalently  >j/>0( /s)  , where  0 is  the  usual  order 
symbol.  Recall  that  the  outer  expansion  (Eq.  1.1-5)  fails  for 
ii<0(z).  Thus,  the  inner  and  outer  expansions  share  an  over- 
lapping region  of  convergence  0(z)  < v < 0(/e).  This  overlap 
implies  that  a uniformly  valid  "composite"  expansion  can  be 
constructed  from  the  inner  and  outer  expansions. 
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In  general  a composite  expansion  can  be  constructed 
from  an  n-term  outer  expansion  and  m-term  inner  expansion, 
where  n and  m are  positive  integers.  For  simplicity  1-term 
expansions  will  be  used.  First,  the  common  part  (cp)  of  the 
inner  and  outer  expansions  is  determined  by  writing  the  inner 
expansion  in  outer  variables  and  expanding  for  small  s. 


2 1/2 

R ) ' 


r-1 ) 2 - u2]1/2 


(1.1-3) 


The  same  result  is  obtained  by  writing  the  1-term  outer  ex- 
pansion in  inner  variables  and  expanding  for  small  c.  Thus, 
the  inner  and  outer  expansions  "match"  in  the  overlapping 
region  of  convergence.  A uniformly  valid  composite  1-term 
expansion  can  ue  obtained  by  either  the  multiplicative  or 
additive  rules  (Ref.  4). 


Tc(r,i|0  = TiTQ/cp 


T (r,-a)  = T±  + TQ  - cp 


(1.1-9) 


(1.1-10) 


The  additive  rule  appears  more  frequently  in  textbooks  (Refs. 
2,  3),  but  the  multiplicative  rule  is  more  convenient  for 
this  example.  The  multiplicative  rule  gives 


( r2-t-l-2rcoso/ ) 2 


2 2 1/2 

[(i+zr+R-] 1/2 


(l.i-n) 


2 1 , 9 

v,  r - j.-2rcos;i/ ) 


(z2  + R2)1/2 

r 2 9 1 , 0 

C(1  ±z)  V R-]  - 


(1.1-12) 


At  low  altitudes,  the  outer  expansion  and  common  part  cancel 
each  other  in  Eq.  1.1-11  so  that  the  composite  solution  is 
dominated  by  the  inner  term.  At  high  altitudes,  the  inner 
expansion  and  common  part  cancel  each  other  in  Eq.  1.1-12 
and  the  outer  term  dominates. 
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Errors  associated  with  the  various  expansions  for  the 
potential  -field  of  the  point  mass  are  illustrated  in  Figure 
1.1-2.  The  depth  of  the  point  mass  is  319  km,  corresponding 
to  e = 0.05.  The  mass  of  the  disturbing  body  is  adjusted  so 
that  an  undulation  of  1 meter  is  generated  directly  over  the 
point  mass.  Figure  1.1-2  shows  that  the  error  in  the  1-term 
inner  expansion  (flat-earth  approximation)  is  zero  directly 
over  the  point  mass  (^=0),  rises  to  1 cm  at  p=4  degrees,  and 
rises  to  1 cm  again  at  b=180  degrees.  The  2-term  inner  expan- 
sion is  more  accurate  than  the  1-term  inner  expansion  in  the 
vicinity  of  the  point  mass  ('Ik30  degrees)  but  suffers  the  same 
error  growth  at  large  angles.  Conversely,  the  1-term  outer 
expansion  is  relatively  accurate  at  large  angles  ('p>30  degrees) 
but  inaccurate  near  the  point  mass.  Thus,  the  inner  and  outer 
expansions  are  complementary:  they  provide  accuracy  in  the 
"near-field"  and  "far-field",  respectively. 


3 30  SO  30  120  150  180 

ANGLE  p (DEGREES 


Errors  in  Various  Expansions  of  the  Undulation 
Generated  by  a Shallow  Point  Mass  (c=0.05,  D= 

319  km).  The  Undulation  Directly  Over  the  Point 
Mass  is  1 Meter. 
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Figure  1.1-2 


A frequency  interpretation  of 
pansions  can  also  be  made,  as  follows, 
corresponding  to  Eq . 1.1-5  is 


the  inner  and  outer  ex- 
The  Legendre  series 


T0(r’*>  = l “571  caPn(cos-4;) 
n=o  r 

cQ  = (l-e)a  - u(l-ne+ ) 


(1.1-13) 

(1.1-14) 


The  expansion  ( Eq . 1.1-14)  is  accurate  for  low-order  coefficients, 
but  fails  for  n>l/s.  In  other  words,  the  outer  expansion  (Eq. 
1.1-5)  provides  only  long-wavelength  information.  The  inner  ex- 
pansion (Eq.  1.1-7)  has  the  opposite  property:  it  provides  only 
short-wavelength  information.  Together  the  two  expansions  des- 
cribe the  gravity  field  accurately  over  all  wavelengths. 


For  comparison,  the  errors  in  the  various  expansions 
are  illustrated  in  Figure  1.1-3  for  a deeper  point  mass  (.£=0.20, 
D=1270  km).  As  before,  the  mass  of  the  disturbing  body  is  ad- 
justed so  that  an  undulation  of  1 meter  is  generated  directly 
over  the  point  mass.  The  error  in  the  1-term  inner  expansion 
(flat-earth  approximation)  rises  to  4 cm  at  b = l 8 degrees  and 
again  at  ^=180  degrees.  The  2-term  inner  expansion  and  1-term 
outer  expansion  provide  improved  accuracy  in  the  near  field 
and  far  field,  as  already  seen  in  Figure  1.1-2. 

At  this  point  it  is  worthwhile  to  review  the  key  ideas. 
Five  formulas  have  been  given  for  the  potential  field  of  a point 
mass:  the  exact  formula  (Eq.  1.1-1)  and  four  expansions  - 

Legendre  (Eq.  1.1-4),  outer  (Eq.  1.1-5),  inner  (Eq.  1.1-7),  and 
composite  (Eq.  1.1-11).  The  round-earth  formulas  (Eqs.  1.1-1, 
1.1-4  and  1.1-5)  are  useful  at  high  altitudes  (r>l+s)  but  in- 
accurate at  low  altitudes.  (For  example,  when  D=10  km  and 
altitude  = 5 km  the  exact  formula  and  single-precision  (7-place) 
arithmetic  yield  only  1-place  accuracy.)  The  inner  expansion 
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*•31082 


Figure  1.1-3  Errors  in  Various  Expansions  for  the 

Undulation  Generated  by  a Shallow  Point 
Mass  (£=0.20,  D=1270  km).  The  undulation 
Directly  Over  the  Point  Mass  is  1 Meter. 


(Eq.  1.1-7)  is  useful  at  low  altitudes  but  fails  to  converge 
at  large  angles  (ij>>0(/s)).  The  composite  expansion  (Eq.  1.1-11) 
enjoys  the  best  properties  of  round-  and  flat-earth  geodesy: 
it  is  accurate  everywhere  (all  altitudes  and  angles).  Gener- 
alizing from  this  example,  the  following  lesson  can  be  drawn. 
Round-  and  flat-earth  formulas  are  useful  for  representing 
long-  and  short-wavelength  gravity  disturbances,  respectively. 
Composite  expansions  provide  a uniformly  valid  representation 
over  all  wavelengths. 


It  should  be  emphasized  that  the  usual  flat-earth 
approximation  is  often  adequate.  The  accuracy  improvements 
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afforded  by  inner  and  composite  expansions,  such  as  Eqs . 1.1-7 
and  1.1-11,  are  typically  rather  small  (less  than  one  percent 
in  the  example,  Figure  1.1-2).  However,  the  expansions  are 
still  quite  useful  because  they  provide  an  error  analysis  of 
the  flat-earth  approximation  and  a systematic  method  for  im- 
proving the  approximation,  if  such  improvements  are  needed. 

The  Fourier  transform  of  the  1-term  inner  expansion 
is  defined  as 

00 

?,(<*»,,  2)  - j'f  T,  ( x , y , z ) el('Jlx^27)dxdy  (1.1-15) 

Inasmuch  as  the  potential  field  of  the  point  mass  is  isotropic, 
the  two-dimensional  Fourier  transform  (Eq.  1.1-15)  can  be 
written  as  a one-dimensional  Hankel  transform  (Ref.  5). 


00 


T1(oj,z) 

= 2t  / RT,  (R,  z ) J (ujR)dR 
o 1 ° 

00 

(1.1-16) 

J ('jR)dR  (1.1-17) 

* 4 / - ) > r 9 211/7 

0 [(l+z)  + IT  J17- 

2 

2 2 

0) 

3 U1  + w2 

(1.1-18) 

This  integral 

is  given  by  Ref.  6,  thus 

T .(a, z) 

2ttu  -(1+z)oj 

(1.1-19) 

Notice  that  the  Fourier  transform  (Eq.  1.1-19)  of  the  point 
mass  field  is  unbounded  at  zero  frequency. 
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STATISTICAL  EXAMPLE  - ATTENUATED  WHITE  NOISE  ACF 


The  following  model  was  recently  proposed  for  the 
gravity  disturbance  potential  autocovariance  function  (acf) 
(Ref.  7)  : 


M r , i ) 

£ 3 ( 2— £ ) 2 a2  [r2  — ( 1— £ ) 4 

( -|  0 

[l-(l-£)4' 

||^r2+(I-£)4  - 2(l-£)2  r cos  tf> 

3/2  ^ 

where 

r 

- r1ra/S* 

e - D/Re 

(1.2 

Here  r^.r^,,  and  are  the  radii  and  angular  separation  of  the 

points  being  correlated.  The  model  contains  two  free  para- 
o 

meters  (a  ,D)  that  are  fitted  to  the  data.  The  Legendre  ex- 
pansion of  Eq . 1.2-1  is 


t ( r , •!) ) 


-:3(2-£l3 

Fl-Cl-s)4’  (l-i 


00 


l (2n+l) 
n=o 


il±l 

r 


2 


a+1 

? ( cosb ) 
n 

(1.2-3) 


For  modeling  local  (short  wavelength)  gravity  distur- 
bances, the  characteristic  distance  D is  small  relative  to 
the  radius  of  the  earth.  (Typically,  D=50  km  and  s=0.008). 
Therefore,  Eqs . 1.2-1  and  1.2-3  are  ill-conditioned  and  slowly 
converging,  respectively.  Expanding  in  powers  of  the  parameter 
e yields 


2s2  a2(r2-l) 

?(r,'40 5 372 

0 (r  +l-2r  cos  ; V ; ' 


lxc  2 ( 5r2-l ) - 2 ( l+3r2 ) r cos  v 
(r2-l)(r2+l  - 2r  cos  w) 


(1.2.4) 
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This  (outer)  expansion  diverges  when  r<l+0(s)  and  #<0(s).  As 
in  the  earlier  example,  inner  variables  (z,R)  are  defined  by 
Eq.  1.1-6.  The  inner  expansion  is 


d»(a,z)  - 


4C2y-z)  a‘ 


2 + 3/2 


’ 2+z ) +•  R 


o o oo 

1+e  z~(2-i-zP  + 2 ( 3-z"  ) R" 

2 (2+z)  ( (2+z)2  + R2] 


(1.2-5) 


As  before,  the  first  term  of  the  inner  expansion  is  the  famil- 
iar flat-earth  approximation.  The  inner  expansion  fails  at 
high  altitudes  (z>2 / z)  and  large  shift  distances  #>0(1). 

A composite  expansion  will  now  be  formed  from  the  1- 
term  inner  and  outer  expansions.  The  common  part  (cp)  is 


o 

4 Q 2 


CP  ’ (z2*a2)3/2 


2 0 

4s2o~(r-l) 

r,  . .2  . , 21  3/2 
(r-1)  to  ' 


(1.2-6) 


and  the  multiplicative  composite  expansion  is 


, . Rr-l)2  + #2]  3/2  r+1  4(  2+z )c2  / 1 o 

$ r ( * > # ) ''"  o 3/2  0 r 0 013/9 

c (r  +l-2r  cos#)3'2  ,2  ( 2+z ) -+R“  j 3/ - 


2 0 0 
2s2Q-‘(r  -1) 


2-z  (zW)3/2 


( r2+l-2r  cos# ) 3 ^ 2 z 


2 + ,2-13/2 


(1.2-8) 


>+z)-  + R 


The  Fourier  transform  of  the  1-term  inner  expansion 
is  given  by  Eqs . 1.1-15,  1.1-16  and  1.2-5,  thus 


$i(+,z)  = 3“(  2+z ) 'i2  / - --  J0  (-iR)dR 

0 i_(  2-z  ) " - R2  j 3/  - 


Cl- 2-9) 


L 


I 


This  integral  is  given  in  Ref.  6,  hence 


. , . 0 2 -( 2+z  )cu 

^(uj.z)  = 8va  e 


(1.2-10) 


Notice  that  unlike  the  previous  example  (Eq.  1.1-19),  the 
Fourier  transform  (Eq.  1.2-10)  is  bounded  at  zero  frequency. 
These  transforms  (Eqs.  1.1-19  and  1.2-10)  are  referred  to  in 
later  sections  of  the  paper. 

For  simplicity,  both  the  deterministic  and  statis- 
tical examples  are  based  on  isotropic  gravity  fields.  For 
nonisotropic  fields,  the  appropriate  transformation  between 
outer  variables  (r,il;,a)  and  inner  variables  (x,y,z)  is 


x =*  (<|>/s)  sin  a 
y = (b/e)  cos  a 
z = (r-l)/s 


(1.2-11) 

(1.2-12) 

(1.2-13) 


The  origin  of  the  inner  coordinate  system  is  on  the  surface  of 
the  sphere  at  the  center  of  the  local  field  (e.g.,  directly 
above  the  point  mass).  The  spherical  angles  'i>  and  a denote 
the  angular  separation  and  azimuth  of  a point  relative  to  this 
origin  (Figure  1.2-1). 
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2.  UPWARD  CONTINUATION  VIA  MATCHED  ASYMPTOTIC  EXPANSIONS 


Poisson's  integral  for  upward  continuation  above  a 
spherical  earth  is  well  known  (Ref.  8). 


2 i 2tt  7T  f f 

T(r,9,X)  = T (9  ,X  ) sine  d9  dX 

4 IT  o o 


(2-1) 


2 2 

l - r + 1 - 2r  cosW 


cosW  53  cos9cos9  +•  sin9sin9  cos  (X  - X) 


(2-2) 

(2-3) 


If  many  evaluations  of  this  integral  are  needed,  it  is  more 
efficient  to  use  transforms: 

oo 

T(  r , 9 , X)  = l -j-rr  Y ( 9 , X ) (2-4) 

n— o r x 

2iT  TT  r 

Yq(  9 , X ) = / / T°(9',x')  Pn(cosW)  sin9  ' d9  ' dX  ' 

o o 

(2-5) 


However,  the  harmonic  expansion  (Eq.  2-4)  converges  too  slowly 
for  modeling  local  ( short -wavelength)  gravity  disturbances. 
Local  disturbances  can  be  extrapolated  aloft  with  the  flat- 
earth  version  of  Poisson's  integral. 


T(x,y 


00 


T°(u,v)  dudv 

Q x-u ) 2 + ( y-v ) 2 +z2]  3/2 


(2-6) 


Flat-earth  upward  continuation  can  also  be  done  using  Fourier 
transforms  (Appendix  A). 
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T(X,  7, 2)  = // 

4tt“  — » 


/ ; - ( y ! » 2 ^ e 


-^2  -i(a1x+-a,7)rfii  , , 


T°(ai,tu2)  = //  T°(x,y)  el(uJlx+aJ2y)  dxdy 


(2-7) 

(2-S) 


However,  the  flat-earth  formulas  (Eqs.  2-6  to  2-8)  are 
only  valid  in  the  vicinity  of  the  origin  (small  x,  y and  z) . 
Therefore,  neither  the  round-  nor  flat-earth  formulas  provide 
an  accurate,  unif ormly-valid  method  for  upward  continuation 
of  local  gravity  disturbances. 

In  this  section,  an  inner  expansion  is  derived  for 
upward  continuation  of  local  gravity  data.  The  nonuniformity 
of  this  expansion  at  large  distances  is  corrected  by  matching 
with  an  outer  expansion.  The  resulting  composite  expansion 
is  uniformly  valid  yet  retains  the  computational  advantages 
of  the  flat-earth  formulas  (Eqs.  2-6  to  2-8).  In  other 
words,  the  composite  expansion  makes  it  possible  to  use  Fourier 
transforms  for  upward  continuation  and  obtain  results  that 
are  accurate  throughout  the  entire  space  (r,9,\)  outside  the 
sphere.  For  simplicity,  the  derivation  is  carried  out  for  the 
special  case  of  an  isotropic  gravity  field.  The  nonisotropic 
expansions  are  inferred. 


The  disturbance  potential  satisfies  a Dirichlet  oroblem 


outside  a snhere. 


o 1 

7 Lr  + 


cop  T 3 n | Kr<  = 

2 b I — 30  < :j)  < » 

r ' 


(2-9) 


T = Tu(*;e) 


(2-10) 


T°(<|>;e)  = T°(  ’^+2t  ; e ) 


(2-11) 


Subscripts  denote  differentiation  (3T/3r=Tr).  It  is  assumed 
that  the  boundary  potential  T°  depends  on  a small  parameter 
e,  and  can  be  expanded  in  a power  series 


T (<|/;e)  - T°(R)  * sT°(R)  + 


(2-12) 


where  R=^/t.  An  example  of  such  a power  series  is  the  shallow 
point  mass,  where 


T (ili’f)  - ±L= i + i 

C*'  5 ( 1+R2 ) 1//2  2 


(2-13) 


according  to  Eq.  1.1-7.  In  more  practical  situations  such  as 
upward  continuation  of  local  gravity  data,  the  dimension  of 
the  local  field  is  D,  t=D/R  , and  Eq . 2-12  truncates  to  a 

_ V 

single  term  (T°a0,  n=2, 3 ,...).  (Such  an  example  is  given  in 
the  next  section). 


An  asymptotic  inner  expansion  is  sought  of  the  form 


T1(r,’i»;e)  - T^z.R)  + sT2(z,R)  + ... 


(2-14) 


where  z=(  r-1 ) /s . ?rhen  Eq.  2-14  is  substituted  into  Eq.  2-9, 
the  following  sequence  of  boundary-value  problems  is  ob- 
tained . 


First  order  problem 


’ + i T 

1RR  R 1R 


0<Z<» 


.2-15) 


21 


Tx  = (R) 


(2-16) 


Second-order  oroblem 


: + i T,  + T2  = -22  T - 2T 

2RR  R 2R  22  X2Z  Z 


T2  = T°(R) 


(2-17) 


(2-18) 


The  solutions  are  found  by  Hankel  transforms  (Ref.  5) 


First  order  solution  (isotropic  case) 


T^(z,R)  = 2t  j 


1 (•*  ( jj  ) e Jg  ( Ro ) da) 


(u)  = 2t  f RT°(R)  Jq(Ro)  dR 

'o 


(2-19) 


(2-20) 


Second  order  solution  (isotroaic  case) 

T2(z,R)  = - | (Tx  - z T1  ) 

z 


j * (a»)  e*,JJZJ0(Ro)  doj 


(2-21) 


T?(w)  » 2 rr  | RT°(R)  J (Ro)  dR 

2 jQ  2 o 


(2-22) 


;f  the  boundary  condition  is  nonisotropic,  the  Hank el  trans- 
forms are  replaced  by  Fourier  transforms. 
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First  order  solution  (nonisotropic  case) 


T t ( x , 7 , z ) = -K  / / f?  (w  uOe'^e 


—oiZ  -i(  j.  ) . 

1FuJ2)e  e 1 2y 'd&j.,  dtu2 


(2-23) 


/ / T°(x,7)  dxdy  (2-24) 


Second  order  solution  (nonisotropic  case) 


T2(x,7,z)  = - f (T1+zT1  ) 

21 


- _ -i(  'jj-  x+u0y  ) 

+ — ^ //  T2  («  «2)e-w  e du1du2 

47T  — « 

(2-25) 


T^Cwi.Wg)  = 


.*  o K^x+u^y) 

JJ  T„(x,y)  e “ “ dxdy 


(2-26) 


A comparison  of  Eqs.  2-7  and  2-23  reveals  that  the  first  term 
of  the  inner  expansion  (Eq.  2-14)  is  the  familiar  flat-earth 
approximation.  This  confirms  that  Eqs.  1.2-11  through  1.2-13 
is  an  appropriate  choice  of  inner  variables. 

The  Fourier  transform  of  the  two-term  inner  expansion 


Ti(oj1,(u2>z)  - TjCa^.u^)  e z 1 - (1-cuZ) 


(2-27) 


g,  rpO  / n _ — 2 

~ ^2  ,JJl,a)2  e 


In  many  cases,  such  as  real  data  processing,  the  boundary  con- 
dition is  only  one  term: 
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T°(’p,i;s)  = T°(x,y) 


T;(x,y)  = 0 


n=2 , 3 , 


(2-28 


( 2-29) 


This  is  a special  case  which  simplifies  the  solutions  (Eas. 
2-19  to  2-27).  In  particular,  the  parameter  t disappears 
from  the  solution,  when  the  solution  is  written  in  dimensional 
form.  The  dimensional  variables  are  primed: 


Q 

X 

II 

X 

y = yD  z = zD 

(2-30) 

1 

= w./D 

» ? ’ 2 ' 2 
,JJ2  = w * ^-j  ^2*" 

(2-31) 

olution  (Eq 

. 2-27)  becomes 

r . » > 

Ti(Jjl’(U2  ’ 2 

* r 

V ^ t!°  ( , , ' \ “'-U  Z 

, 1 

1 - — — ( 1 -n  z ) 1 

L 2He  J 

(2-32) 

Equation  2-32  provides  a very  simple  recipe  for  upward  contin- 
uation of  potential  field  data.  First,  the  data  is  Fourier  trans- 
formed to  obtain  the  two-dimensional  array  T°('jj,  ,j9).  This 
array  is  multiplied  by  the  appropriate  factors  in  Eq.  2-32 
and  inverse  transformed  to  obtain  the  final  results.  In  the 
limit  of  large  earth  radius.  Eq . 2-32  degenerates  into  the 
usual  flat-earth  formula.  Of  course,  the  inner  expansion 
(Eq.  2-14)  is  only  valid  near  the  origin  (small  x,y,z).  An 
outer  expansion  is  needed  to  remove  this  limitation. 


The  foregoing  analysis  provides  a method  for  compu- 
ting the  inner  expansion  using  Hankel  or  Fourier  transforms. 
Direct  integrals  can  also  be  used,  as  shown  in  Appendix  B. 
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However,  the  Fourier  transforms  are  more  efficient  than 
numerical  integration. 


Ltii 


The  outer  expansion  will  now  be  constructed.  On  an 
outer  scale,  the  local  boundary  condition  Eq.  2-10  shrinks  to  a 
point  as  the  parameter  s goes  to  zero.  Therefore,  the  outer 
expansion  consists  of  singularity  functions  that  satisfy  two 
requirements:  they  must  obey  Laplace's  equation  and  match  the 

inner  expansion.  Examples  of  singularity  functions  that 
satisfy  Laplace's  equatioa  have  already  been  studied  in  the 
previous  sections,  namely 


T(r,’l)  = u(r2  + l-2r  cosid)*1^ 


(2-33) 


and 


T(r,v) 


u(r  -1) 


(r2  + l-2r  cost;)3/2 


(2-34) 


To  match  the  inner  and  outer  expansions,  the  1-term 
inner  expansion  is  written  in  outer  variables  and  expanded 
for  small  £ . 


cp 


- 


J 


1 <«)  e~*(r_1)/sJ oUa/e)  da  (2-35) 


A change  of  variable  (u^u/e)  yields 

2 >oo 


cn 


z 
YT 


1 0 


u T^(eu)  e"(r*I)u  J (*u)  du  (2-36) 


The  long-wavelength  properties  of  the  local  field  are  paramount 
as  £ goes  to  zero.  If  the  transform  is  bounded  at  zero 
frequency  then  Eq . 2-36  becomes 


co  = 


2t 


• 00 

1(o)J0  u e“(r"1)uJo(’iiu)  du 


(2-37) 
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3ere  T^(z,R)  denotes  the  first  term  of  the  inner  expansion, 
namely  Eq . 2-19. 

Returning  now  to  Eq . 2-36,  suppose  that  the  transform 
is  unbounded  at  zero  frequency  as  in  the  first  example 
(Eq.  1.1-19).  If  the  local  field  satisfies 


lim  '.a  T^(uj)  = A (2—42) 

O 

where  A is  a constant,  then  Eq.  2-36  becomes 


Ia  this  case,  the  appropriate  singularity  for  the  outer  expansion 
is  Eq.  2-33  with 


U 3 A £ / ( 2 it) 


(2-45) 


The  composite  expansion  is 


t (,.♦>  - -4sd>.Viva 


(r2+l-2r  costy)^2 


T1(z,a) 


(2-46) 


Inasmuch  as  the  local  boundary  condition  Eq,  2-10  shrinks  to  a 
point  as  the  parameter  £ goes  to  zero,  the  outer  expansion 
and  common  parr  are  the  same  for  both  isotropic  and  noniso- 
tropic  local  fields.  Therefore,  the  appropriate  composite 
expansion  for  nonisotropic  local  fields  is 


Tc(r,i|))a) 


C(r-l)Wl3/2 

(r2+i-2r  cosb)3/2 


r+1 

2 


T1(x,y,z) 


(2-47) 


if  the  following  integral  is  bounded  and  non-zero: 

30 

//  T°(x,y)  dxdy 


(2-48) 


However  if  the  integral  (Eq.  2-48)  is  unbounded  and 


lim 

0j-*-O 


w //  T"(x,y)  e1^  1 


i(u-x+u„y) 


dxdy 


('2-49) 


where  A is  a constant,  then  the  composite  expansion  is 


Tc(r,^,a)  ~ 


[(r-1)2  + 'i;2]1/2  „ , 

— 2 V2  T1(x,y,z) 

(r  +l-2r  cost!/)”' 


(2-50) 


27 


FT 


The  reader  should  keep  in  mind  that  the  composite  expansions, 
such  as  Eq.  2-41,  are  uniformly  valid.  Thus,  upward  contin- 
uations of  local  fields  can  be  done  using  flat-earth  formulas 
and  when  the  flat-earth  results  are  multiplied  by  the  appro- 
priate factor  in  Eq.  2-41,  they  become  valid  for  all  altitudes 
and  angles  (all  r,ip). 

When  processing  data,  the  local  field  will  usually 
have  bounded  non-zero  energy  at  zero  frequency,  i.e.,  Eq . 2-48 
is  bounded  and  non-zero.  Therefore,  Eq . 2-47  is  the  appropriate 
composite  expansion.  A notable  exception  to  this  rule  is  a 
local  field  due  to  one  or  more  point  masses.  In  this  case, 

Eq . 2-43  is  unbounded  and  Eq . 2-50  is  the  proper  expansion. 

Both  the  disturbance  potential  and  the  disturbance 
potential  acf  satisfy  a Dirichlet  problem  outside  a sphere 
(Ref.  9).  Therefore,  expansions  for  upward  continuation  of 
acfs  can  be  inferred  from  the  previous  results.  For  upward 
continuation  of  the  disturbance  potential,  (r,  v , a)  denote 
normalized  radius  (r^r'/R  ),  angle,  and  azimuth.  For 
upward  continuation  of  acfs,  (b, a)  denote  the  angle  and 
azimuth  between  the  points  being  correlated  and 

r = r1r2/Rg  (2-51) 

where  r^  and  r9  are  the  radii  of  the  two  points. 
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3. 


LIMITED  DATA  EXAMPLE 


The  earlier  examples  are  instructive  because  they 
provide  a simple  exposition  of  the  key  notions: 


• For  purposes  of  modeling  local  gravity 
disturbances,  matched  expansions  and 
Fourier  transforms  are  more  accurate  and 
efficient  than  spherical-earth  formulas 
and  Legendre  transforms. 

• The  first  term  of  the  inner  expansion 

is  the  familiar  flat-earth  approximation. 
Improved  accuracy  is  obtained  by  comput- 
ing additional  terms. 

• The  nonuniformity  of  the  inner  expansion 
at  large  distances  is  corrected  by  outer 
and  composite  expansions. 


Also,  the  section  on  upward  continuation  shows  that  the  inner 
expansion  can  be  computed  with  Fourier  transforms. 


Unfortunately,  neither  of  the  earlier  examples 
address  the  most  common  type  of  local  gravity  field,  namely, 
the  gravity  field  over  a limited  patch  of  the  earth's  surface. 
A prototype  problem  is  to  compute  the  exterior  potential  from 
data  within  the  patch.  The  boundary  condition  is 


T ( , a ; s ) 


t(x,y) 

0 


-1  < x,y  < 1 
otherwise 


(3-1) 


Presumably,  data  outside  this  region  is  continued  upward  by  seme 
other  means.  For  example,  data  outside  the  region  might  be 
continued  upward  using  spherical-earth  (Legendre  transform) 
methods.  Alternatively,  if  the  earth's  surface  is  divided  into 
convenient  blocks  of  data,  each  block  can  be  continued  upward 
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using  Four*  methods.  Thus,  it  is  sufficient  to  consider  the 
problem  of  a single  block.  To  simplify  the  following  example, 
a circular  area  is  used,  and  the  boundary  potential  is  con- 
stant over  the  area. 


T°(4>;s) 


t ii<e 

0 0»>£ 


(3-2) 


The  more  realistic  case  of  a square  patch  and  nonuniform 
potential  can  be  treated  by  essentially  the  same  methods.  In 
particular,  the  1-term  inner  and  composite  expansions  for  the 
general  case  are  Eqs . 2-23  and  2-47.  Inasmuch  as  the  boundary 
condition  is  isotropic,  the  exact  solution  Eq . 2-1  can  be 
reduced  to  a single  integral  (Ref.  9) 


T(r^)  , (r2-l)  rff  T° (’!>;£ ) sine  E(v)  d0 

^ 0 (a-b)  / a+b 


(3-3) 


where 


a 


r i-l-2r 


cosbcosS 


(3-4) 


b = 

E(  v) 

2 


2r  sia’^sinG 

rff/2 


/l-v2  sia2x  dx 


2b/ ( a+b ) 


(3-5) 

(3-6) 


(3-7) 


Here  E 
itmd . 


denotes  the  complete  elliptic  integral  of 
Substituting  Eq . 3-2  into  3-3  yields 

T(r  rJ)  a Izf-U.  - r 3inJ_EivJ_^e 

■0  (a-b)  va-^o 


the  second 


(3-3) 
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The  inner  expansion  will  now  be  constructed, 
boundary  condition  Eq . 3-2  becomes 


T°  (R) 


I15 

»o 


(R)  = 0 

n 


R<1 

R>1 

n-2,3, . . 


The 


(3-9) 


(3-10) 


The  first-order  inner  solution  is  Eq . 2-19  where 

r1 


T^(u») 


2 if  t 


i R J (ujR)  dR 
0 

This  integral  is  given  by  (Ref.  6): 

T^(oj)  = 2fT(  t/u ) J1(uj) 

Therefore,  the  first-order  inner  solution  is 

f00 

T,  ( z , R)  =*  t ! e”,'oZ  J (;uR)  j.(u)  du 
i0  ° 1 

This  integral  cannot  be  expressed  in  terms  of  elementary 
functions.  The  second-order  inner  solution  is 


(3-11) 


(3-12) 


(3-13) 


T2(z.R)  - — | (Tx  + z Tl  ) 


(3-14) 


Inasmuch  as  the  transform  T^  is  bounded  at  zero  frequency,  the 
1-term  outer  expansion  is  Eq . 2-34  with 


c ^ — o 2 

U - f*  T^(0)  « £2t/4 


The  1-term  composite  expansion  isEq.  2-41. 


(3-15) 


For  the  special  case  p=0 , both  the  exact  solution 
(Eq.  3-3)  and  matched  expansions  can  be  evaluated  in  terms  of 
elementary  functions.  The  exact  solution  (Eq.  3-8)  becomes 

T(r,o)  = [*1  - (r-1 ) / (r2+l-2  r cose)I/2j  (3-1 C) 

The  first-order  inner  solution  can  be  evaluated  from  either 
the  Hankel  transform  Eq . 3-13  or  the  direct  integral  Eq . B-l : 


T-,  (z,o)  = t (1  - z/Zl+O 


The  1-term  outer  expansion  Eq . 2-34  for  '1=0  is 


T0  (r.o)  - <e2t/4> 


and  the  1-term  multiplicative  composite  expansion  is 


(3-17) 


(3-13) 


(3-19) 


The  second-order  inner  solution  is 
T2  ( z ,0)  - -t(z/2)  | i - z 


2+z‘ 


ci^h3/2 


(3-20) 


Notice  that  for  a local  gravity  disturbance  (s<<1), 
the  exact  solutions  (Eqs.  3-8  and  3-16)  are  numerically  ill- 
conditioned.  The  spherical-earth  formulas  do  not  provide  an 
accurate,  efficient  method  for  upward  continuation.  However, 
the  matched  expansion  (Eq.  3-19)  is  accurate,  efficient,  and 
uniformly  valid.  Furthermore,  the  matched  expansion  can  be 
easily  improved,  if  necessary,  by  adding  an  additional  term 
to  the  inner  and/or  outer  expansion.  If  more  accuracy  is 
needed  in  the  near  field,  another  term  is  added  to  the  inner 
expansion.  Similarly,  another  term  in  the  outer  expansion 
improves  accuracy  in  the  far 


field. 


4. 


STOKES ' INTEGRAL  VIA  MATCHED  ASYMPTOTIC  EXPANSIONS 


Stokes'  integral  relates  gravity  anomalies  to  the 
disturbance  potential  on  the  sphere 

R 2rr  tt 

T( 9 , X ) 38  -r”  / / Ag(ti,ct)  S(iO  sinpdpdf.  (4-1) 
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where  S is  Stokes'  function  (Ref.  8).  If  many  evaluations 
of  this  integral  are  needed,  a spherical  harmonic  transform 
is  useful . 


T( 9 , X ) 


l 

n=*o 


ifn 

n-1 


(4-2) 


» 2u+l  / jT  dg(a]\')  PQ(cosb)  sin 9 d9  dA  (4-3) 

^ * ' o o 

However,  the  harmonic  series  (Eq.  4-2)  converges  too  slowly  for 
modeling  local  (short-wavelength)  gravity  anomalies.  The 
disturbance  potential  due  to  local  gravity  anomalies  can  be 
computed  using  the  flat-earth  version  of  Stokes’  integral 

(Ref.  10). 


T(x, y) 


_D 

2tt 


Ag( u , v)  du  dv 
((x-u)2  + (y-v)2]1^2 


(4-4) 


The  flat-earth  computation  can  also  be  carried  out  using  courier 
transforms  (Appendix  A). 


T(x,y) 


j f (l/oj)Ag  ('-^  > ^2  ^ 


9-i(Ulx+«l27)dUi4oj2 

(4-5) 


Ag(ui1,  m2) 


00 

/ / Ag(x.y)  e 

-.00 


i(Ulx+w2y)  dxdy 


(4-6) 
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However,  the  flat-earth  formulas  (Eqs.  4-4  through  4-6)  are  only 
valid  in  the  vicinity  of  the  origin  (small  x,y).  Therefore, 
neither  the  round-  nor  flat-earth  formulas  provide  an  accurate, 
unif ormly-valid  method  for  computing  the  disturbance  potential 
from  local  gravity  anomalies. 

In  this  section,  inner,  outer,  and  composite  expansions 
are  derived  for  computing  the  disturbance  potential  from  local 
gravity  anomalies.  As  in  the  previous  section,  the  derivation 
is  carried  out  for  the  isotropic  case  and  corresponding  non- 
isotrcpic  results  are  inferred. 


Stokes'  integral  is  the  solution  to  the  differential 


equation 


rT  + 2T  = -r(r, 1) 
r 


(4-7) 


The  function  f (f=rrdg)  satisfies  a Dirichlet  problem  outside 
the  sDhere 


‘rr  r 


1 j,  , cotu  - _ n 

.2  ' r2  * 'jj  " 0 


(4-3) 


f - f°  (*;£)-  R 4g 


(4-9) 


f°(il;;s)  3 f°  ('ii+2ir ; s ) 


(4-10) 


As  in  a previous  section  asymptotic  expansions  are  assumed 
of  the  form  Eq . 2-14  and 


f(r,’4;s)  ~ f^z.R)  f9(z,R)  +• . . . 


(4-11) 


- f°(R)  -e  z°2  (?.)+.  . . 


( il-12  ) 
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The  following  sequence  of  boundary-value  problems  is  obtained. 
First-order  problem 


T * -ef 
1,  1 
z 


(4-13) 


Second-order  problem 


T2z  = " £ f2  - 2T1Z  - 2T1 


(4-14) 


The  solutions  are  readily  found,  inasmuch  as  the  Dirichlet 
problem  was  treated  in  a previous  section. 

First-order  solution  (isotropic  case) 


Tl(2,R)  - 


(m)  e~^z  Jo(0)R)  dm 


(4-15) 


f^(m)  - 2,rr  j Rf°  (R)  J0(ojR)  cLR 


(4-16) 


Second-order  solution  (isotropic  case) 


T2(z,R)  - - zTx  + -g-y-  j 


(1/m)  f^  (ai)  e”^"  Jq(uR)  dm 


+ 4z  f f,(w)  e'*2  J (rnR)  dm 

2TJ0  2 ° 


(4-17) 


f^(m)  - 2.TT  j R f 2 ( R)  JQ(wR)  dR 


(4-18) 


At  zero  altitude  (z»0,  r » 1)  we  have 
TX(R)  - f *1  («)  J0(uR)  dm 


(4-19) 


35 


A comparison  of  Eqs . 4-5  and  4-21  reveals  that  the  first  term 
of  the  inner  expansion  (Eq.  2-14)  is  the  familiar  flat-earth 
approximation.  This  confirms  that  Eqs.  1.2-11  through  1.2-13 
is  an  appropriate  choice  of  inner  variables. 

The  outer  expansion  will  now  be  constructed  bv  a 
matching  procedur.  as  was  done  in  a previous  section.  The 


1-term  inner  expansion  (Eq.  4-15)  is  written  in  outer  variables 
and  expanded  for  small  t . 


CP 


C - TO/,.')  a 

2x  ' o ^ 


-ij(r-ly/s.  j ( -jjiu / s ) du 


(4-25) 


A change  of  variable  (u=u/c)  yields 

.2 


cp 


r £ *?«»>  - 


-(r-l)u  jo^u)  du 


(4-26) 


Assuming  that  the  integral 


f^(0)  = 2 t R f°  (R) 


dH 


(4-27) 


is  bounded  and,  non-zero,  Eq . 4-26  becomes 

cp  = 2*7  £^(0)/  e"(r'1)u  J0(H»u)  du 
o 


-2_o  f 2 9 / 2 

frj-  f^(0)/  ! ( r-1 ) 2 + / Jx/ 2 


( 4-2S ) 


(4-29) 


The  implied  singularity  in  the  outer  expansion  is  Eq.  2-33 
where 


U * e2  f^(0)/(  2tt) 


The  composite  expansion  is 


(r  „ _ Ur-1)2 * 1^1 

:(r,i0  / _2  , , xl/2 


(r  + l-2r  cosv)' 


T1(z,R) 


(4-30) 


(4-31) 


• ] 

1 

■; 


At  zero  altitude  (r=l,  z-Q)  we  have 


T„(1/)  - 


/2( 1-cosy; ) 


1/2 


(H) 


(4-32) 
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a(?/2)  csc( y/2 ) T ( R) 


(4-33) 


The  corresponding-  nonis'otropic  result  ’is 


T ( li , ct ) = ( ii/2 ) • csc(  i»  / 2 ) T.,  (xj) 


(4-34) 


, The  term  ( 'b/2  )csc (i/2 ) in  Eq . 4-34  is  a correction  factor 
applied  to  the  flat-earth  results.  This  factor  is  near  unity 
— for  example,  when  ^=10  and  60  degrees,  the  factor  is  1.001 
and  1.047,  respectively.  Equations  4-33  and  4-19  constitute  a 
unif ormly-valid  solution  to  Stokes'  problem.  Local  gravity 
anomalies  are  processed  using  the  flat-earth  formulas;  when 
the  results  are  multiplied  by  the  appropriate  factor  in  Eq . 
4-33,  they  become  valid  for  all  angles  (u). 

Returning  now  to  Eq . 4-26,  suppose  that  the  integral 
(Eq.  4-27)  is  zero.  Then  Eq . 4-26  becomes 


CO  — z 


~ J ( j u ) du 
c 


(4-35) 


s 3 3 ( r-1) 
((r-1)2  +'l;2]  3/2 


(4-36) 


where 


-lim  . 2 -jo 


(R)  J, (uR)  dR 


(4-37) 


The  implied  singularity  in  the  outer  expansion  is 
Eq.  2-34  with 

u - (-:3/2)  3 


(4-3S) 
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The  composite  expansion  is 


T„(r.uO 


9 9 7/9 

(r-ir  + ■!>-]  3/~  r-1 


2 


(r-fl-2r  cosip) 


77^72  2 


T1(z,R) 


At  zero  altitude  (r»l,  z=0)  we  have 


Tc(<|»)  - ( TX> 3 / 8 ) csc3(^/2)  T1(R) 


The  corresponding  nonisotropic  result  is 


Tc(i|»,a)  ~ (ip3/3)  csc3(ip/2)  T^x.y) 


Note  that  Eqs . 4-34  and  4-41  apply  when  the  integral 
(Eq.  4-27)  is  non-zero  and  zero,  respectively- 


(4-39) 


(4-40) 


(4-41) 
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VENING-MEINESZ  INTEGRAL  VIA 
MATCHED  ASYMPTOTIC  EXPANSIONS 


The  Vening-Meinesz  integral  relates  gravity  anomalies 
to  vertical  deflections  (Ref.  8). 


in(9,X)>  4irC  ~S 


O O 


( cosa ) dS 

a)  l \ -j—  s m’pdw da 

(sina) 


(5-1) 


For  dealing  with  local  ( short-wavelength ) gravity  anomalies,  tb 
flat-earth  version  is  convenient  (Ref.  10). 


f?(x, 7) 

\ o ( X , 7 ) 


1 

2"G  J ! 


dg(u,  v)  (y~v) 

5 b-  • < \dudv  (5-2) 

C(x-u)2  + (y-v)2]3/2  |x-u| 


The  flat-earth  computations  can  also  be  carried  out  with  Fourier 
transforms 


5(x, y) 
n ( x , y ) 


4tt  G ~=° 


TT  [ f J W e'i(ill^27)^2 


( 5-3 ) 


where  Ag  is  given  by  Eq . 4-6.  However,  the  flat-earth  formulas 
(Eqs.  5-2  and  5-3)  are  only  valid  in  the  vicinity  of  the 
origin  (small  x,y). 


Asymptotic  expansions  that  enjoy  the  best  properties 
of  the  round-  and  flat-earth  formulas  can  be  derived  by  dif- 
ferentiation of  the  results  in  the  previous  section.  In 
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IS 


4tt2GD 


i / 


(■-1/'a»)  ( uj1  , ;jj2 ) 9'l("Jlx4“'iJ27)du1du; 


The  associated  formulas  for 


are  obvious. 


(3-10) 


The  outer  expansion  will  now  be  constructed  by  the  same 
matching  procedure  used  in  the  previous  sections.  The  1-term 
inner  expansion  (Eq . 5-7)  is  written  in  outer  variables  and 
expanded  for  small  s . 


CD 


-( x/R)e 
2tGD 


/a f^(oj)  e”'JJ^r~1)/,~  J 1((*np/s ) du 


A change  of  variable 


( u=»'jj / t ) yields 


(5-11) 


cp 


-( x/R)e> 
2fGD 


j uf  t (eu)  e 
o 


-(r-l)u 


J1(iu) 


du 


(5-12) 


Assuming  that  Eq . 4-27  is  bounded  and  non-zero 

(x/R)e3tO  i 

cp  = s o.at;  1 - 1 ( 0 ) 5 b 3 72 

2tGD  1 [ ( r-1)2  + b2]  ' 


(5-13) 


The  implied  singularity  in  the  outer  expansion  is 


j , _ ,t.  \ ur  sin  b 


(r  -^l-2r  cosiu ) 


3/2 


(5-14) 


where 


u = i^£Vt(0) 

2ttGD  1 


( 5-15) 
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The  1-term  multiplicative  composite  expansion  for  both  the 
isotropic  and  nonisotropic  cases  is 

(x,y)  - ( ^ 2 / 8 ) sinU;  csc3(t|;/2)<5  (x,y)  (5-16) 

c X1 

Returning  now  to  Eq.  5-12,  suppose  that  the 
integral  (Eq.  4-27)  is  zero.  Then  Eq.  5-12  becomes 


cp  3 - B f u2e"(r~1)uJ,  ('^u)  du 

GD  ^ 1 


(5-17) 


(x/R)s4  B — 3(r~1)’£ 

GD  |^(  r-1 ) 2 4-  -j;2j5/2 


(5-18) 


The  implied  singularity  for  the  outer  expansion  is 

5 ( r , i|0  - -Gursint;  (r2-l) 
x (r2+l-2r  cos^)3/2 


(5-19) 


wher‘ 


(x/R)s' 

GD 


B 


(£-20) 


The  1-term  multiplicative  composite  expansion  for  both  the 
isotropic  and  nonisotropic  cases  is 

5 (x.7)  - (t(>4/32)  sinb  csc5(-J;/2)  5 (x,y)  (5-21) 

4 c X1 

The  foregoing  expansions  provide  the  vertical  deflec- 
tion components  (5^,5  ) in  the  local  x,y  coordinate  system. 
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A transformat ion  is  necessary  to  convert  these  components 
into  north  and  east  deflections  ( C , n ) • This  transformation 
depends  on  the  latitude  and  longitude  of  the  origin  of  the 
local  coordinates,  as  well  as  the  position  (x,y)  of  the  point 
where  the  deflection  is  sought.  Although  tedious,  such  trans- 
formations are  elementary,  and  are  not  presented  here. 
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6. 


DISCRETE  FOURIER  TRANSFORMS 


The  foregoing  theory  is  suitable  for  analytical  model 
such  as  the  earlier  examples.  For  numerical  models,  the  dis- 
crete Fourier  transform  is  more  useful . The  disturbance 
potential  T and  its  transform  T are  both  treated  as  discrete 
and  periodic.  The  Fourier  transform  .(Eq.  2-3)  becomes 


M/2-1  N / 2-1 

t°  =.  y > if .» 

an  H--11/2  1--N/2 


2Tri(mk/M+ni/N) 


(6-1) 


where  the  integers  k,2,,m,n  are  defined  by 

x » k ix  y = i iy 


(6-2) 


oj1  = mft. 


m2  m ai22 


(6-3) 


n. 


2tt 


a. 


2* 


(6-4) 


“1  - MAX ' 2 NAy 

For  simplicity,  M and  N are  assumed  to  be  even  integers.  The 
inverse  transform  for  upward  continuation  (Eq.  2-7)  becomes 

M/2-1  N/2-1 

T . _L  V Y t°  a-2Tri(mk/M  + q*/n>  (R 

kl  Nil  L 2-  ran  w (.o--) 

m*-M/2  n~-N/2 


where 


2 2 
*1  + J2 


2„2  2„2 

mAn.  + n ^2 


(6-6) 


Fast  Fourier  transform  algorithms  are  normally  based  on  one- 
sided transforms  whereas  Eqs . 6-1  and  6-6  are  two-sided.  The 


appropriate  one-sided  transforms  can 

the  summation  operators,  inasmuch  as 

T,  . and  transform  T are  oeriodic  . 
K-t  mn 


be  obtained  by  shifting 
both  the  potential  field 
The  one-sided  transforms 


M-l  N-l 

T°  = y y T°  P 2Tri(mk/M  + nl/X) 
-mn  ,4-  LkZ  e 


(6-7) 


M-l  X-l 


mkZ 


y y 

NM 


m=o  n=o 


T° 

mn 


-wz  a -2fi(mk/M  + ni/X) 


where 


2 2 o 2 2 

'jj  = p + a ^2 


p = | m - M/2  -M/2 


q = In  - N/2 ! -N/2 


(6-8) 


(6-9) 

(6-10) 

(6-11) 


As  an  example  of  upward  continuation  using  discrete 
Fourier  transforms,  consider  the  potential  field  of  a point 
mass , 


D * depth  of  point  mass  * 319  km  (6-12) 

lx  » Ay  * D/2  = 160  km  (6-13) 

z » 319  km  (6-14) 

Case  1:  M = X = 32  (6-15) 

Case  2:  M = X = 64  (5-16) 


As  in  the  earlier  example,  the  mass  of  the  disturbing  body  is 
selected  so  that  an  undulation  of  1 meter  is  obtained  directly 
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over  the  point  mass.  The  point  mass  is  located  directly  below 
the  center  of  the  square  grid.  In  each  case,  the  potential 
field  of  the  point  mass  is  computed  at  zero  altitude.  These 
gridded  values  are  continued  aloft  using  the  discrete  Fourier 
transforms  (Eqs.  6-7  and  6-8).  The  results 'are  compared  with  the 
exact  flat-earth  solution  to  assess  the  errors  introduced  by 
the  discrete  transforms.  Profiles  of  the  errors  are  shown  in 
Fig.  6-1.  For  the  smaller  grid,  the  error  is  2 to  3 cm,  where- 
as the  larger  grid  reduces  the  errors  to  0.2  to  0.3  cm.  The 
larger  grid  is  obviously  a better  approximation  of  the  infinite- 
length  transforms  (Eqs.  2-7  and  2-8). 

Errors  introduced  by  the  flat -earth  approximation  and 
discrete  Fourier  transforms  can  be  compared  by  examining  Figs. 
1.1-2  and  6-1.  Figure  1.1-2  shows  that  the  flat-earth  approxima- 
tion introduces  an  error  of  T cm.  Thus,  the  error  introduced 
by  discrete  transforms  is  larger  than  the  flat-earth  error  for 
the  smaller  grid  (M=N=32)  and  smaller  than  the  flat-earth 
error  for  the  larger  grid  (M=N=64).  These  results  are  not 
surprising  because  the  dimensions  of  the  smaller  and  larger 
grids  are  5100  x 5100  km  and  10,200  x 10,200  km,  respectively, 
whereas  the  earth's  mean  radius  is  6371  km. 

It  is  interesting  to  compare  the  frequency  content 
of  the  errors  introduced  by  the  flat-earth  approximation  and 
discrete  transforms.  The  flat-earth  error  (Fig.  1.1-2)  is 
basically  binodal:  a short-wavelength  spike  appears  at 
’1=4  degrees  and  a long-wavelength  error  appears  at  '1=180 
degrees.  The  discrete  transform  error  is  also  bimodal: 
short-  and  long-wavelength  errors  arise  from  the  discrete 
spacing  and  finite  length  of  the  data.  Of  course,  these 
short-wavelength  errors  decay  quickly  as  altitude  increases. 
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DISTANCE  FROM  CENTER  OF  GRID  (km) 


Figure  6-1 


Errors  in  discrete  Fourier  Transform 
Upward  Continuation  of  the  Undulation 
Generated  by  a Shallow  Point  Mass. 

The  Undulation  Directly  Over  the 
Point  Mass  is  1 Meter. 


The  point  mass  example  is  convenient  for  comparing 
errors  introduced  by  the  flat-earth  approximation  and  finite- 
length  transforms.  However,  this  example  is  somewhat  benign 
because  the  potential  field  tapers  smoothly  to  zero  at  the 
edges  of  the  finite  grid.  Real  gravity  fields  do  not  taper 
to  zero  at  the  edges  and  discontinuities  appear  in  a periodic 
extension.  In  such  cases,  the  errors  introduced  by  finite- 
length  transforms  are  more  serious  than  portrayed  in  the  point 
mass  example  (Fig.  6-1  )•  However  techniques  are  available  for 
mitigating  such  "edge  effects"  (Refs.  11,  14). 


For  readers  who  are  unfamiliar  with  the  fast  Fourier 
transform,  some  remarks  about  computer  time  are  appropriate. 
The  computer  time  needed  for  this  algorithm  tends  to  be  quite 
minimal,  and  increases  at  a relatively  modest  rate  (MN  log  MN) 
as  the  dimension  (MxN)  of  the  grid  increases  (Ref.  12). 


48 


The  foregoing  example  problem  required  5 seconds  of  computer 
time  on  the  IBM  370/158  digital  computer,  which  included  four 
Fourier  transforms  (two  for  M=N=32  and  two  for  M=N=64). 


7. 


CONCLUSIONS 


Flat -earth  formulas  provide  efficient  descriptions 
of  local  ( short -wavelength)  gravity  disturbances.  The  accuracy 
of  this  approximation  can  be  improved  by  the  method  of  matched 
asymptotic  expansions.  The  inner  and  outer  expansions  remove 
the  short-  and  long-wavelength  errors  of  the  flat-earth 
formulas,  respectively.  The  resulting  composite  expansions 
provide  uniformly  valid  descriptions  of  all  wavelengths  over 
all  altitudes  and  distances.  The  additional  computations 
needed  to  achieve  these  accuracy  improvements  are  typically 
insignificant  because  the  inner  expansion  is  computed  using 
Fourier  transforms  and  the  outer  expansion  consists  of  simple 
singularity  functions  such  as  Eqs . 2-33  and  2-34. 

If  discrete  rather  than  continuous  Fourier  transforms 
are  used,  additional  errors  arise  due  to  the  finite  spacing 
and  length  of  the  grid.  However,  these  errors  can  be  reduced 
by  using  a finer  grid  spacing,  larger  grid,  and  "edge  effect" 
compensation  techniques. 
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APPENDIX  A 

FOURIER  TRANSFORMS  OF  POISSON  AND  STOKES ' INTEGRALS 


The  flat-earth  Poissoa  Integral  (Eq.  2-6)  can  be  viewed 
as  a linear  system  with  the  boundary  potential  T°  and  exterior 
potential  T as  the  input  and  output.  The  impulse  and  fre- 
quency response  functions  are 


h(x,y) 


2ir(x2+y2+z2)3/2 


(A-l) 


* 2ff 


fl  ,„222,-3/2 


i(w,  x + u0y ) 


i 

.00 


(x4'+yi+2“)  “ s 1 2 dxdy  (A-2) 


Inasmuch  as  the  impulse  response  is  isotropic,  the  Fourier 
transform  (Eq.  A-2)  can  be  written  as  a Hankel  transform 


H(u>1,ui2) 


where 


R(R2  -k  z2)"3/2  Jq(u»R)  dR 


2 _ 2 , 2 


(A-3) 


(U  * v dig 


R2  = x2  + y2 


( A-4 ) 
(A-5) 


Equation  A-3  is  the  same  transform  that  appears  in  the 
statistical  example  (Eq.  1.2-9)  hence 


H(  ui) 


— ‘jjZ 


(A-6) 


which  yields  Eq.  2-7. 
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= 


The  flat-earth  Stokes'  integral  (Eq.  4-4)  can  be  analy- 
zed similarly.  The  impulse  and  frequency  response  functions  are 


fcL(x.y)  = 


2f  (x2^2)1^2 


(A-7) 


i o i/o  x™f.u  ) 

( x“+y  ; "1/  ^ e 41  dxdv 


(A-8) 


= D 


J (uR)  dH 
o 


(A-9) 


This  is  the  same  transform  that  appears  in  the  deterministic 
example  (Eq..  1.1-17)  hence 

H(w)  - D/o)  ( A-10 ) 


which  yields  Eq . 4-5. 
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APPENDIX  3 

DIRECT  INTEGRAL  FORM  OF  INNER  EXPANSIONS 


In  the  text,  the  inner  expansions  such  as  Eq . 2-14  are 
expressed  in  terms  of  Hankei  and  Fourier  transforms,  for  example, 
Eqs . 2-19,  2-20  and  2-23,  2-24.  These  solutions  can  also  be  ex- 
pressed as  direct  integrals.  For  example  Eqs.  2-19  and  2-20  are 
replaced  by  (Ref.  13). 


T1(2,a) 


2z  7 

TT  ^ 

o 


T°(u)  E(p)  udu 
C ( R-u) 2+z2]  C(R-u)2-z2] 1/2 


(B-l) 


wt\ere  E is  the  complete  elliptic  integral  of  the  second  kind 
(Eq.  3-6)  and 


„2  4Ru 

0 = 6 2 

(R+u)“+z^ 


(B-2 ) 


Equation  3-1  effectively  replaces  two  integrals  by  one.  Sim- 
ilarly, the  second-order  inner  solution  (Eqs.  2-21  and  2-22) 
can  be  expressed  as 


T2(z,R) 


- (z/2) 


(T1+zT1  ) 

2 

T2(u)E(p)  udu 
[ (R-u)2+z2]  ((R^u)2+z2] 1/2 


( B-3 ) 


For  the  nonisotropic  case,  Eqs.  2-23  and  2-24  are  replaced  by 
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(x,  7. 2)  = 5T- 


T^(u,v)  dudv 


; ; 

—co 


9 9 21/ 9 

[( x-u )"+ (y-v)“  -2'"  ] ' “ 


(3-4) 


Similarly,  Eqs . 2-25  and  2-26  are  replaced  by 


T0 (x , y , z ) = - 


(z/2)  (T1+zT1  ) 
z 

® To ( u , v ) dudv 

r t 2 

J l 


z 

2ir 


C(x-7)2-^(7-v)  .z2)3/2 


(3-5) 


In  the  Stokes'  solution,  Eqs.  4-21  through  4-24  are 
replaced  by 


Tl(x,y)  - 27  // 


f-^Cu.v)  dudv 
-«  ( ( x-u)2+( y-v) 2] 1/2 


(3-6) 


and 


T2(X'7)  = 47  fl 


f^(u,v)  In  [ ( x-u) 2+( y-v ) 2 ] dudv 


» f g( u < v)  dudv 

+ 27  -»  [(x-u)2+(y-v)2] 1/2 


(B-7) 


In  order  for  the  first  integral  in  Eq . B-7  to  be  bounded,  it 
is  required  that 


if  f°(x,y)  dxd7  = 0 (3-8) 

—CO  -*■ 

The  disadvantage  of  these  direct  integrals  is  that 
numerical  integrations  are  needed  for  each  term  of  the  inner 
expansion.  When  Fourier  transforms  are  used,  numerical  inte- 
grations are  avoided  altogether. 
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